🏗️ Constructing single-cell trajectories with Monocle3🧜

1 🔃 Load required libraries

library(ggplot2)
library(htmlwidgets)
library(igraph)
library(monocle3)
library(plotly)

set.seed(12345)

2 🧬 Constructing single-cell trajectories

During development, in response to stimuli, and throughout life, cells transition from one functional “state” to another. Cells in different states express different sets of genes, producing a dynamic repetoire of proteins and metabolites that carry out their work. As cells move between states, they undergo a process of transcriptional re-configuration, with some genes being silenced and others newly activated. These transient states are often hard to characterize because purifying cells in between more stable endpoint states can be difficult or impossible. Single-cell RNA-Seq can enable you to see these states without the need for purification. However, to do so, we must determine where each cell is in the range of possible states. For a comparison of different trajectory inference methods, see Saelens et al. (2019) (Saelens et al. 2019).

Monocle 3 (Cao et al. 2019) is the latest version of the Monocle toolkit. Monocle introduced the strategy of using RNA-Seq for single-cell trajectory analysis (Trapnell et al. 2014). Rather than purifying cells into discrete states experimentally, Monocle uses an algorithm to learn the sequence of gene expression changes each cell must go through as part of a dynamic biological process. Once it has learned the overall “trajectory” of gene expression changes, Monocle can place each cell at its proper position in the trajectory. You can then use Monocle’s differential analysis toolkit to find genes regulated over the course of the trajectory, as described in the section Finding genes that change as a function of pseudotime. If there are multiple outcomes for the process, Monocle will reconstruct a “branched” trajectory. These branches correspond to cellular “decisions”, and Monocle provides powerful tools for identifying the genes affected by them and involved in making them. You can see how to analyze branches in the section Analyzing branches in single-cell trajectories.

The workflow for reconstructing trajectories is very similar to the workflow for clustering, but it has a few additional steps. To illustrate the workflow, we will use another C. elegans data set, this one from Packer & Zhu et al. (Packer et al. 2019). Their study includes a time series analysis of whole developing embyros. We will examine a small subset of the data which includes most of the neurons.

3 📥 Loading the data

We will load the data as we did with the L2 data. This dataset contains an expression matrix, cell metadata, and gene annotations.

expression_matrix <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_expression.rds"))
cell_metadata <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_colData.rds"))
gene_annotation <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_rowData.rds"))

cds <- new_cell_data_set(expression_matrix,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)

4 ⚙️ Pre-process the data

Pre-processing works exactly as in clustering analysis. This time, we will use a different strategy for batch correction, which includes what Packer & Zhu et al did in their original analysis.

Note: Your data will not have the loading batch information demonstrated here, you will correct batch using your own batch information.

cds <- preprocess_cds(cds, num_dim = 50)
cds <- align_cds(cds, alignment_group = "batch", residual_model_formula_str = "~ bg.300.loading + bg.400.loading + bg.500.1.loading + bg.500.2.loading + bg.r17.loading + bg.b01.loading + bg.b02.loading")

Note that in addition to using the alignment_group argument to align_cds(), which aligns groups of cells (i.e. batches), we are also using residual_model_formula_str. This argument is for subtracting continuous effects. You can use this to control for things like the fraction of mitochondrial reads in each cell, which is sometimes used as a QC metric for each cell. In this experiment (as in many scRNA-seq experiments), some cells spontanously lyse, releasing their mRNAs into the cell suspension immediately prior to loading into the single-cell library prep. This “supernatant RNA” contaminates each cells’ transcriptome profile to a certain extent. Fortunately, it is fairly straightforward to estimate the level of background contamination in each batch of cells and subtract it, which is what Packer et al did in the original study. Each of the columns bg.300.loading, bg.400.loading, corresponds to a background signal that a cell might be contaminated with. Passing these colums as terms in the residual_model_formula_str tells align_cds() to subtract these signals prior to dimensionality reduction, clustering, and trajectory inference. Note that you can call align_cds() with alignment_group, residual_model_formula, or both.

5 🎨 Reduce dimensionality and visualize the results

Next, we reduce the dimensionality of the data. However, unlike clustering, which works well with both UMAP and t-SNE, here we strongly urge you to use UMAP, the default method:

cds <- reduce_dimension(cds)
p1 <- plot_cells(cds, label_groups_by_cluster=FALSE,  color_cells_by = "cell.type")
ggsave("figures/monocle3_trajectory_umap.pdf", p1, width=12, height=10)
ggsave("figures/monocle3_trajectory_umap.jpg", p1, width=12, height=10)

Figure 1: UMAP of C. elegans embryonic cells colored by cell type

As you can see, despite the fact that we are only looking at a small slice of this dataset, Monocle reconstructs a trajectory with numerous branches. Overlaying the manual annotations on the UMAP reveals that these branches are principally occupied by one cell type.

As with clustering analysis, you can use plot_cells() to visualize how individual genes vary along the trajectory. Let’s look at some genes with interesting patterns of expression in ciliated neurons:

ciliated_genes <- c("che-1",
                    "hlh-17",
                    "nhr-6",
                    "dmd-6",
                    "ceh-36",
                    "ham-1")

p2 <- plot_cells(cds,
           genes=ciliated_genes,
           label_cell_groups=FALSE,
           show_trajectory_graph=FALSE)
ggsave("figures/monocle3_trajectory_ciliated_genes.pdf", p2, width=12, height=10)
ggsave("figures/monocle3_trajectory_ciliated_genes.jpg", p2, width=12, height=10)

Figure 2: Expression of ciliated neuron genes along the trajectory

We will learn how to identify the genes that are restricted to each outcome of the trajectory later on in the section Finding genes that change as a function of pseudotime.

6 🔬 Cluster your cells

Although cells may continuously transition from one state to the next with no discrete boundary between them, Monocle does not assume that all cells in the dataset descend from a common transcriptional “ancestor”. In many experiments, there might in fact be multiple distinct trajectories. For example, in a tissue responding to an infection, tissue resident immune cells and stromal cells will have very different initial transcriptomes, and will respond to infection quite differently, so they should be a part of the same trajectory.

Monocle is able to learn when cells should be placed in the same trajectory as opposed to separate trajectories through its clustering procedure. Recall that we run cluster_cells(), each cell is assigned not only to a cluster but also to a partition. When you are learning trajectories, each partition will eventually become a separate trajectory. We run cluster_cells() as before.

cds <- cluster_cells(cds)
p3 <- plot_cells(cds, color_cells_by = "partition")
ggsave("figures/monocle3_trajectory_partitions.pdf", p3, width=12, height=10)
ggsave("figures/monocle3_trajectory_partitions.jpg", p3, width=12, height=10)

Figure 3: UMAP of C. elegans embryonic cells colored by partition

7 🕸️ Learn the trajectory graph

Next, we will fit a principal graph within each partition using the learn_graph() function:

cds <- learn_graph(cds)

p4 <- plot_cells(cds,
           color_cells_by = "cell.type",
           label_groups_by_cluster=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE)

ggsave("figures/monocle3_trajectory_graph.pdf", p4, width=12, height=10)
ggsave("figures/monocle3_trajectory_graph.jpg", p4, width=12, height=10)

Figure 4: UMAP of C. elegans embryonic cells with learned trajectory graph

This graph will be used in many downstream steps, such as branch analysis and differential expression.

8 ⏳ Order the cells in pseudotime

Once we’ve learned a graph, we are ready to order the cells according to their progress through the developmental program. Monocle measures this progress in pseudotime.

What is pseudotime?

Pseudotime is a measure of how much progress an individual cell has made through a process such as cell differentiation.

In many biological processes, cells do not progress in perfect synchrony. In single-cell expression studies of processes such as cell differentiation, captured cells might be widely distributed in terms of progress. That is, in a population of cells captured at exactly the same time, some cells might be far along, while others might not yet even have begun the process. This asynchrony creates major problems when you want to understand the sequence of regulatory changes that occur as cells transition from one state to the next. Tracking the expression across cells captured at the same time produces a very compressed sense of a gene’s kinetics, and the apparent variability of that gene’s expression will be very high.

By ordering each cell according to its progress along a learned trajectory, Monocle alleviates the problems that arise due to asynchrony. Instead of tracking changes in expression as a function of time, Monocle tracks changes as a function of progress along the trajectory, which we term “pseudotime”. Pseudotime is an abstract unit of progress: it’s simply the distance between a cell and the start of the trajectory, measured along the shortest path. The trajectory’s total length is defined in terms of the total amount of transcriptional change that a cell undergoes as it moves from the starting state to the end state.

In order to place the cells in order, we need to tell Monocle where the “beginning” of the biological process is. We do so by choosing regions of the graph that we mark as “roots” of the trajectory. In time series experiments, this can usually be accomplished by finding spots in the UMAP space that are occupied by cells from early time points:

p5 <- plot_cells(cds,
           color_cells_by = "embryo.time.bin",
           label_cell_groups=FALSE,
           label_leaves=TRUE,
           label_branch_points=TRUE,
           graph_label_size=1.5)
ggsave("figures/monocle3_trajectory_time_bin.pdf", p5, width=12, height=10)
ggsave("figures/monocle3_trajectory_time_bin.jpg", p5, width=12, height=10)

Figure 5: UMAP of C. elegans embryonic cells colored by embryo time bin

The black lines show the structure of the graph. Note that the graph is not fully connected: cells in different partitions are in distinct components of the graph. The circles with numbers in them denote special points within the graph. Each leaf, denoted by light gray circles, corresponds to a different outcome (i.e. cell fate) of the trajectory. Black circles indicate branch nodes, in which cells can travel to one of several outcomes. You can control whether or not these are shown in the plot with the label_leaves and label_branch_points arguments to plot_cells. Please note that numbers within the circles are provided for reference purposes only.

Now that we have a sense of where the early cells fall, we can call order_cells(), which will calculate where each cell falls in pseudotime. In order to do so order_cells() needs you to specify the root nodes of the trajectory graph. If you don’t provide them as an argument, it will launch a graphical user interface for selecting one or more root nodes.

# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, time_bin="130-170"){
  cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin)

  closest_vertex <-
  cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <-
  igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
  (which.max(table(closest_vertex[cell_ids,]))))]

  root_pr_nodes
}

cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))

In the above example, we just chose one location, but you could pick as many as you want. Plotting the cells and coloring them by pseudotime shows how they were ordered:

p6 <- plot_cells(cds,
           
           color_cells_by = "pseudotime",
           label_cell_groups=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE,
           graph_label_size=1.5)

ggsave("figures/monocle3_trajectory_pseudotime_interactive.pdf", p6, width=12, height=10)
ggsave("figures/monocle3_trajectory_pseudotime_interactive.jpg", p6, width=12, height=10)

Figure 6: UMAP of C. elegans embryonic cells colored by pseudotime

Note that some of the cells are gray. This means they have infinite pseudotime, because they were not reachable from the root nodes that were picked. In general, any cell on a partition that lacks a root node will be assigned an infinite pseudotime. In general, you should choose at least one root per partition.

It’s often desirable to specify the root of the trajectory programmatically, rather than manually picking it. The function below does so by first grouping the cells according to which trajectory graph node they are nearest to. Then, it calculates what fraction of the cells at each node come from the earliest time point. Then it picks the node that is most heavily occupied by early cells and returns that as the root.

# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, time_bin="130-170"){
  cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin)

  closest_vertex <-
  cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <-
  igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
  (which.max(table(closest_vertex[cell_ids,]))))]

  root_pr_nodes
}
cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))

Passing the programatically selected root node to order_cells() via the root_pr_nodes argument yields:

p7 <- plot_cells(cds,
           color_cells_by = "pseudotime",
           label_cell_groups=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE,
           graph_label_size=1.5)

ggsave("figures/monocle3_trajectory_pseudotime_programmatic.pdf", p7, width=12, height=10)
ggsave("figures/monocle3_trajectory_pseudotime_programmatic.jpg", p7, width=12, height=10)

Figure 7: UMAP of C. elegans embryonic cells colored by pseudotime

Note that we could easily do this on a per-partition basis by first grouping the cells by partition using the partitions() function. This would result in all cells being assigned a finite pseudotime.

9 🌿 Subset cells by branch

It is often useful to subset cells based on their branch in the trajectory. The function choose_graph_segments allows you to do so interactively.

cds_sub <- choose_graph_segments(cds)

10 🧊 Working with 3D trajectories

You can also perform trajectory analysis in 3D.

cds_3d <- reduce_dimension(cds, max_components = 3)
cds_3d <- cluster_cells(cds_3d)
cds_3d <- learn_graph(cds_3d)

cds_3d <- order_cells(cds_3d, root_pr_nodes=get_earliest_principal_node(cds))

# Get the UMAP coordinates and partition information
umap_coords <- reducedDims(cds_3d)$UMAP
cell_partitions <- partitions(cds_3d)

# Create a data frame for plotting
plot_df <- data.frame(
  UMAP_1 = umap_coords[,1],
  UMAP_2 = umap_coords[,2],
  UMAP_3 = umap_coords[,3],
  Partition = factor(cell_partitions)
)

# Create and display interactive 3D plot with plotly
# Create the 3D plot
fig <- plot_ly(plot_df, 
              x = ~UMAP_1, 
              y = ~UMAP_2, 
              z = ~UMAP_3,
              color = ~Partition,
              type = "scatter3d",
              mode = "markers",
              marker = list(size = 3)) %>%
  layout(scene = list(
    xaxis = list(title = "UMAP_1"),
    yaxis = list(title = "UMAP_2"),
    zaxis = list(title = "UMAP_3")
  ))

# Display the plot directly
fig

References

Cao, J., M. Spielmann, X. Qiu, X. Huang, D. M. Ibrahim, A. J. Hill, F. Zhang, et al. 2019. “The Single-Cell Transcriptional Landscape of Mammalian Organogenesis.” Journal Article. Nature 566 (7745): 496–502. https://doi.org/10.1038/s41586-019-0969-x.

Packer, J. S., Q. Zhu, C. Huynh, P. Sivaramakrishnan, E. Preston, H. Dueck, D. Stefanik, et al. 2019. “A Lineage-Resolved Molecular Atlas of c. Elegans Embryogenesis at Single-Cell Resolution.” Journal Article. Science 365 (6459). https://doi.org/10.1126/science.aax1971.

Saelens, W., R. Cannoodt, H. Todorov, and Y. Saeys. 2019. “A Comparison of Single-Cell Trajectory Inference Methods.” Journal Article. Nat Biotechnol 37 (5): 547–54. https://doi.org/10.1038/s41587-019-0071-9.

Trapnell, C., D. Cacchiarelli, J. Grimsby, P. Pokharel, S. Li, M. Morse, N. J. Lennon, K. J. Livak, T. S. Mikkelsen, and J. L. Rinn. 2014. “The Dynamics and Regulators of Cell Fate Decisions Are Revealed by Pseudotemporal Ordering of Single Cells.” Journal Article. Nat Biotechnol 32 (4): 381–86. https://doi.org/10.1038/nbt.2859.